# R code for centroids-capitals example
# in CShapes paper

# Nils B. Weidmann <weidmann@icr.gess.ethz.ch>
# December 10, 2008

library(foreign)
library(cshapes)
library(maptools)

# This program requires the Fearon & Laitin (2003) data to be present 
# in the working directory. This file can be obtained from
# http://www.stanford.edu/~jfearon/data/apsr03repdata.zip 

# modify this line, enter your working directory
setwd("/Users/nilsw/Projects/CShapes/R")

# load Fearon&Laitin (2003) dataset
fldata <- read.dta("repdata.dta", convert.factors=F)
fldata$onset[fldata$onset>1] <- 1 # fix coding error
fldata$ccode[fldata$ccode==260 & fldata$year>=1990] <- 255 # update to COW 2008 
fldata <- fldata[fldata$year>=1950,]
fldata$uniqueid <- fldata$ccode*10000 + fldata$year

# gdpenl has many missings. Replace with KSG trade and gdp data.
ksgdata <- read.table("pwt_cow.asc", header=T) 
ksgdata$uniqueid <- ksgdata$statenum*10000 + ksgdata$year
ksgdata <- subset(ksgdata, select=c(uniqueid, gdppc))
fldata <- merge(fldata, ksgdata, by="uniqueid")

fldata <- subset(fldata, 
  !is.na(onset) 
  & !is.na(warl) 
  & !is.na(gdppc) 
  & !is.na(lpopl1) 
  & !is.na(lmtnest) 
  & !is.na(ncontig) 
  & !is.na(Oil)
  & !is.na(nwstate)
  & !is.na(instab)
  & !is.na(polity2l)
  & !is.na(ethfrac)
  & !is.na(relfrac))
  
# compute for each polygon the deviation of the capital from the centroid
cshp <- cshp()
centroids <- coordinates(cshp)
distance <- function(from.x, from.y, to.x, to.y) { 
  from <- matrix(c(from.x, from.y), ncol=2, nrow=1) 
  spDistsN1(from, c(to.x, to.y), longlat = TRUE)
}
cshp$capcentdist <- mapply(distance, centroids[,1], centroids[,2], cshp$CAPLONG, cshp$CAPLAT)

# compute the diagonale used for the relative capital deviation
bboxes <- sapply(seq(1:nrow(as.data.frame(cshp))), function(x) as.numeric(bbox(cshp[x,])))
cshp$diagonale <- mapply(distance, bboxes[1,], bboxes[2,], bboxes[3,], bboxes[4,])

# convert to time series cross-sectional
data.yearly <- cshapes2yearly(cshp, vars=c("capcentdist", "diagonale"), useGW=F)

# merge
data.yearly$uniqueid <- data.yearly$ctrcode*10000 + data.yearly$year
fldata <- merge(fldata, data.yearly, by="uniqueid")
fldata$rel.capcentdist <- fldata$capcentdist / fldata$diagonale

# model 1: absolute deviation alone
m1 <- glm(onset ~ log10(capcentdist), family=binomial, data=fldata)
summary(m1)

# model 2: relative deviation alone
m2 <- glm(onset ~ rel.capcentdist, family=binomial, data=fldata)
summary(m2)

# model 3: All onsets
m3 <- glm(onset ~ log10(capcentdist)+warl+log10(gdppc)+lpopl1+lmtnest+ncontig+Oil+nwstate+instab+polity2l+ethfrac+relfrac, family=binomial, data=fldata)
summary(m3)

# model 4: All onsets
m4 <- glm(onset ~ rel.capcentdist+warl+log10(gdppc)+lpopl1+lmtnest+ncontig+Oil+nwstate+instab+polity2l+ethfrac+relfrac, family=binomial, data=fldata)
summary(m4)

## Plot examples
pdf(file="capcentexamples1.pdf", width=14, height=7)
layout(matrix(1:2, 1, 2, byrow=T))

# DRC
drc <- cshp[cshp$FEATUREID==170,]
plot(drc, lwd=2)
title("Congo (Democratic Republic)", cex=1.5)
points(drc$CAPLONG, drc$CAPLAT, pch=42, cex=5)
points(coordinates(drc), pch=19, cex=2)
lines(c(coordinates(drc)[1], drc$CAPLONG), c(coordinates(drc)[2], drc$CAPLAT))
text(drc$CAPLONG, drc$CAPLAT, "Kinshasa", pos=1, offset=0.9, cex=1.5)
text(coordinates(drc), "Centroid", pos=1, offset=0.9, cex=1.5)


## Nigeria
nigeria <- cshp[cshp$FEATUREID==79,]
plot(nigeria, lwd=2)
title("Nigeria", cex=1.5)
points(nigeria$CAPLONG, nigeria$CAPLAT, pch=42, cex=5)
points(coordinates(nigeria), pch=19, cex=2)
lines(c(coordinates(nigeria)[1], nigeria$CAPLONG), c(coordinates(nigeria)[2], nigeria$CAPLAT))
text(nigeria$CAPLONG, nigeria$CAPLAT, "Abuja", pos=1, offset=0.9, cex=1.5)
text(coordinates(nigeria), "Centroid", pos=3, offset=0.9, cex=1.5)

dev.off()